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SIMULATIONS  OF  GAS  PHASE  DETONATIONS: 
INTRODUCTION  OF  AN  INDUCTION  PARAMETER  MODEL 


I.  Introduction 

The  recent  concentrated  effort  to  understand  those  factors  controlling 
the  detonability  of  gas  mixtures  has  been  triggered  by  the  need  to  assess 
hazards  related  to  both  LNG  spills  and  containment  of  material  in  nuclear 
power  plants.  These  potentially  dangerous  situations  cover  a  range  of 
scenarios  from  unconfined  vapor  cloud  explosions  to  explosions  in  partially 
or  totally  confined  chambers.  The  materials  of  concern  include  mixtures  of 
relatively  light  weight  fuels,  such  as  li^9  CH^,  ^2^2*  etc>  in  either  °2  °r 
mixtures  of  0^  and 

In  this  paper  we  describe  one-  and  two-dimensional  simulations  of  gas 
phase  detonations  of  stoichiometric  mixtures  of  and  CH^  in  air.  Primary 
attention  is  given  to  the  development  and  calibration  of  an  induction  para¬ 
meter  model  to  be  used  in  detonation  simulations.  It  will  be  shown  that 
use  of  such  a  calculation  to  replace  the  integration  of  ordinary  differen¬ 
tial  equations  representing  chemical  kinetics  extends  our  ability  to 
do  calculations  over  physically  meaningful  timescales  in  air  mixtures. 

Simulations  of  gas  phase  detonation  phenomena  have  been  performed  by 
Fickett  and  Wood1  using  one-dimensional  calculations  based  on  the  method 
of  characteristics  in  order  to  study  the  stability  of  the  detonation  front. 
These  tests  were  extended  to  two  dimensions  by  Mader2  who  used  several 

finite  difference  formulations.  All  of  these  calculations,  however, 
Manuscript  submitted  April  17,  1980. 
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assumed  a  one-step  Arrhenius  reaction  with  instantaneous  energy  release. 
Taki  and  Fujiwara3  used  a  two-dimensional  finite  difference  method  and  the 
two-step  reaction  model  used  by  Korobeinikov,  et  al. 4  to  simulate  detona¬ 
tions  in  an  oxy-hydrogen  system.  The  new  contribution  presented  in  this 
paper  is  the  detailed  calibration  of  an  extended  induction  parameter  model 
by  comparison  with  calculations  performed  using  a  detailed  computational 
model  which  includes  a  full  set  of  chemical  reactions.5  After  benchmark 
tests  are  described,  the  induction  parameter  model  is  used  in  one-  and  two- 
dimensional  time-dependent  simulations  of  propagating  and  decaying  detona¬ 
tion  waves. 
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II.  Summary  of  the  Numerical  Models 


The  reactive  shock  models  used  to  perform  the  calculations  described 
below  solve  the  time-dependent  conservation  equations  for  mass,  momentum 
and  energy.6*7  In  addition,  they  may  solve  a  set  of  coupled  ordinary 
differential  equations  describing  the  production  and  loss  of  a  number  of 
chemical  species.  The  convection  terms  are  solved  using  one  variant  of  the 
FCT  algorithm8*9  and  the  ordinary  differential  equations  for  the  chemical 
kinetics  are  solved  using  a  version  of  the  CHEMEQ  algorithm. 10  * 11  An 
extensive  description  of  the  one-dimensional  model5  has  shown  calculations 
of  incident  and  reflected  shocks  in  a  reactive  material.  The  two-dimensional 
model  has  been  used  to  calculate  mixing  and  vortex  formation  at  material 
interfaces12  and  the  effects  of  the  Rayleigh-Taylor  instability  on  the 
surface  of  pellets  being  imploded  by  lasers.13*14 

Integration  of  the  chemical  rate  equations  is  by  far  the  most  expen¬ 
sive  part  of  a  reactive  shock  calculation.  Table  I  summarizes  the  timing 
requirements  for  2000  timesteps  of  a  detonation  simulation  performed  using 
the  one-dimensional  model.  Fifty  chemical  reactions  were  used  to  describe 
the  interactions  among  eight  different  species.  In  the  calculation  of 
detonation  phenomena,  we  are  often  Interested  in  whether  a  particular  igni¬ 
tion  mechanism  with  a  certain  amount  of  initial  input  energy  will  become  a 
self-sustaining  detonation.  In  certain  situations,  such  as  or 

mixtures,  the  characteristic  distances  out  to  which  a  calculation  must  be 
performed  are  a  meter  or  less.  In  such  cases  the  one-dimensional  reactive 
shock  model  may  be  used.  If  multi-dimensional  effects  such  as  cellular 
structure  are  to  be  modelled,  the  two-dimensional  model  may  be  used  at 
great  cost.  However,  for  mixtures  such  as  H^-air  and  CH^-air,  in  which 


the  characteristic  distances  may  be  many  meters ,  the  computational  cost 
would  be  exorbitant* 

Thus  one  new  addition  to  both  the  one-  and  two-dimensional  models  is 
the  introduction  of  a  less  expensive  global  induction  parameter  model 
to  describe  the  chemical  reactions.  The  required  input  for  this  model  is 
a  knowledge  of  the  induction  time  of  a  mixture  as  a  function  of  tempera¬ 
ture  and  pressure,  x(T,P).  Then  a  quantity,  I(x),  is  defined  which  indi¬ 
cates  how  long  the  material  has  remained  at  a  temperature  and  pressure.  This 
quantity  is  convected  with  the  fluid  in  a  spatially  dependent  calculation 
and  it  is  used  to  indicate  when  the  available  chemical  energy  should  be 
released.  The  reaction  is  not  assumed  to  go  immediately  to  completion. 

The  energy  is  released  over  several  timesteps  determined  by  the  ratio  of 
the  sonic  transit  time  to  the  timestep  itself.  The  basic  idea  is  to  keep 
a  whole  computational  cell  from  burning  at  once  to  eliminate  numerical 
"surging”. 

The  adaptive  gridding  method  discussed  by  Oran,  Boris  and  Young5  was 
generalized  for  use  in  the  calculations  presented  below.  To  illustrate 
the  necessity  for  varying  the  spatial  gridding  adaptively,  consider  the 
three  temperature  profiles  shown  in  Fig.  1.  Three  different  resolution 
grids  were  used  to  perform  the  same  physical  calculation,  a  detailed  simu¬ 
lation  of  a  laser-initiated  detonation  of  an  mixture.  The  laser 

was  assumed  to  have  deposited  energy  by  dissociating  or  O2  molecules  in 
a  small  region  at  the  left  hand  side.  The  presence  of  the  radicals  quickly 
triggered  chemical  reactions  and  energy  release.  The  raised  temperature 
and  pressure  then  generated  a  shock  wave  driven  ahead  of  the  flame  front. 


Since  this  is  a  supersonic  problem,  there  is  little  to  be  gained  from 
a  Lagrangian  calculation. The  velocity  of  the  shock  and  detonation 
fronts  will  be  greater  than  the  speed  of  sound  in  the  unreacted  mixture. 
However,  until  the  sustaining  detonation  is  formed,  we  are  in  effect 
dealing  with  a  flame  propagation  problem  behind  a  shock  front.  Thus  we  must 
use  adaptive  gridding  to  obtain  resolution  of  steep  gradients  in  tempera¬ 
ture  and  species  densities.  At  the  short  timescales  of  sonic  transit, 
diffusive  transport  processes  will  do  little  to  smooth  these  gradients. 

In  order  to  cope  with  this  problem  and  achieve  the  good  results  shown 
by  the  solid  curve  in  Fig.  1,  a  relatively  general  adaptive  gridding  method 
was  developed.  As  the  shock  runs  out  ahead  of  the  flame,  one  finely  zoned 
region  of  grid  points  centers  around  the  shock  front  and  an  even  finer  set 
centers  on  the  flame  front.  The  location  of  these  regions  and  the  gradual 
variation  of  cell  size  is  controlled  automatically.  Such  techniques  are  a 
solution  to  the  problem  of  obtaining  the  required  accuracy  while  decreasing 
the  number  of  cells  carried  in  the  calculation.  A  similar  technique  is 
used  in  both  the  one-  and  two-dimensional  models  with  which  the  calculations 
described  below  are  performed. 


III .  The  Induction  Parameter  Model 

To  show  how  the  model  for  t(T,P)  is  derived,  consider  Fig.  2  which 
shows  the  curves  of  induction  time  as  a  function  of  initial  temperature 
and  pressure  for  stoichiometric  mixtures  of  hydrogen  in  air.  This  has 
been  calculated  from  the  chemical  reaction  rate  scheme  constructed  and 
tested  by  Burks  and  Oran.15  In  this  figure  the  induction  time  is  chosen 
to  be  the  time  at  which  the  density  of  the  OH  radical  peaks.  This  is 
one  way  the  induction  time  is  defined  experimentally. 

To  be  most  consistent  with  the  way  in  which  x(T,P)  is  used  in  the 
reactive  flow  model,  the  induction  time  should  be  defined  as  the  time  at 
which  the  energy  of  the  system  begins  to  increase.  This  is  also  rhe  time 
when  the  temperature  begins  to  increase  rapidly.  Figure  3  shows  calcula¬ 
tions  of  temperature  as  a  function  of  time  for  a  stoichiometric  CH^-air 
mixture  at  5  atm  which  is  suddenly  heated  to  temperatures  in  the  range 
from  1800  to  3000 K.  The  reaction  rate  scheme  used  here  was  derived  by 
Gelinas17  and  has  been  shown  to  agree  well  with  measured  induction  times 
for  stoichiometric  mixtures.  Noted  on  these  curves  is  the  time  of  the 
peak  of  the  atomic  oxygen  density. 

Information  such  as  that  summarized  in  Fig.  2  and  shown  in  more 
detail  in  Fig.  3,  or  well  educated  estimates  if  no  detailed  chemical  reac¬ 
tion  rate  scheme  is  available, are  required  as  a  starting  point.  We  note 
that  for  some  of  the  temperatures  and  pressures  characteristic  of  detona¬ 
tions,  there  has  been  no  experimental  data  to  compare  with  the  results  of 
the  chemical  model. 

The  induction  time  curves  are  then  fitted  to  an  equation  of  the  form 
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t(T,P)  -  T.(P)  f  exp  ,  (1) 

where  for  convenience  Pq  *  1  atm.  The  procedure  is  to  first  pick,  one 
curve  at  fixed  pressure  and  guess  a  value  of  T*.  Then  a  least  squares 
fit  is  done  to  determine  the  coefficients,  fm  x*  and  A*.  These  numbers 
are  then  used  to  evaluate  Eq.  (1).  This  process  is  repeated  using  incre¬ 
mented  values  of  T*,  and  the  final  set  of  values  chosen  is  that  set  which 
best  reproduces  the  input  data  for  t.  Calculations  are  made  for  a  range 
of  pressures  from  which  an  analytic  form  for  the  coefficients  is  derived 
as  a  function  of  pressure.  Expressions  for  A*,  T*,  and  m  x*  which  give 
good  fits  to  calculated  induction  times  for  stoichiometric  H^-air  and  CH^- 
air  mixtures  are  given  in  Table  II. 

These  coefficients  reproduce  induction  times  determined  from  chemical 
model  calculations  with  maximum  errors  of  about  20%.  The  largest  errors 
occur  at  low  initial  temperatures  and  at  pressures  greater  than  five 

t 

atmospheres.  We  are  also  required  to  specify  the  energy  released  from 
the  reactions  and  this  information  may  be  determined  from,  for  example, 
the  JANAF  tables.18 

In  the  coefficients  for  both  the  methane-air  and  hydrogen-air  mixtures, 
there  is  a  marked  change  of  shape  in  the  coefficients  at  pressures  greater 
than  a  few  atmospheres.  Our  analysis  of  the  hydrogen-oxygen  chemistry  below 
and  above  this  pressure  shows  that  as  the  pressure  increases,  reactions 
involving  HO2  and  are  more  important  and  HjO  becomes  important  as  a 

reactant.  This  is  consistent  with  the  experimental  observations  of  Voevodsky 
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and  Soloukhin19  and  Meyer  and  Oppenheim20  who  observed  this  behavior  in 
mixtures  in  the  Pressure  and  temperature  range  above  the  second 
explosion  limit  and  below  the  extended  second  limit. 

Finally,  we  note  that  the  form  of  Eq.  (1)  is  similar  to  that  used  by 
both  Meyer  and  Oppenheim20  and  White.21  In  all  three  cases  the  form 
involved  an  inverse  pressure  or  density  dependence  and  an  Arrhenius-type 
exponential.  The  addition  in  Eq.  (1)  of  the  temperature  offset,  T*(P),  in 
the  denominator  of  the  exponential  provides  another  parameter  which  allows 
extra  flexibility  in  describing  the  shape  of  the  induction  time  curves  at 
higher  temperatures.  The  low  temperature  cutoff  itself  does  not  affect  the 
kinetics  in  these  detonation  calculations. 

Korobeinikov  et_  al . 4  and  Taki  and  Fujiwara3  have  given  approximations 
for  the  amount  of  energy  release  and  time  to  equilibration  during  H2~°2 
combustion.  These  approximations  and  our  detailed  chemical  calculations 
show  that,  except  at  high  temperatures,  the  induction  times  are  generally 
long  compared  to  the  characteristic  time  for  energy  release.  The  energy 
release  timescales  are  strongly  influenced  by  the  time  history  of  the 
mixture.  In  detonation  problems  in  which  spatial  and  temporal  variations 
are  important,  both  energy  release  and  induction  time  modelling  are  simi¬ 
larly  limited  in  accuracy.  Higher  order  global  models  could  be  constructed 
which  depend  on  characteristics  of  the  time  variation  of  temperature  and 
pressure,  but  the  inaccuracies  which  occur  when  data  is  forced  to  have  a 
functional  form  outweigh  any  expected  accuracy  gains.  From  the  one¬ 
dimensional  benchmark  comparisons  described  in  the  next  section  we  have 
concluded  that  local  chemical  energy  release  times  can  be  ignored  in  the 
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induction  parameter  model. 

The  simplification  proposed  above  which  does  not  include  an  energy 
release  time  is  further  supported  by  previous  detonation  and  reactive 
shock  calculations.3*5  These  show  that  energy  release  in  a  computational 
cell  is  governed  by  the  time  it  takes  the  deflagration  wave  to  cross  the 
cell  dimensions  rather  than  by  the  energy  release  time  in  the  relatively 
thin  reaction  zone.  In  the  calculations  with  the  induction  parameter 
model  presented  here,  the  available  chemical  energy  is  released  in  one 
acoustic  transit  time  of  the  cell,  which  is  a  good  estimate  of  the  transit 
time  of  the  shock.  Furthermore,  this  approximation  ensures  a  steady 
progression  of  the  computed  reaction  zone  across  the  system.  In  adaptively 
gridded  calculations  this  acoustic  transit  time  algorithm  is  particularly 
important.  The  computational  cells  tend  to  move  with  the  detonation  so  the 
effective  energy  release  rate  is  further  reduced  from  the  peak  chemical 
reaction  rates. 
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IV.  One-Dimensional  Calculations  and  Benchmark  of  the  Induction 

Parameter  Model 

To  benchmark  the  induction  parameter  model  described  above  and  thus 
determine  its  accuracy  and  range  of  validity,  several  tests  were  carried 
out  which  compared  the  results  of  representative  problems  solved  using  the 
full  chemical  reaction  rate  scheme  with  those  solved  using  the  induction 
parameter  model  to  replace  the  detailed  kinetics.  Below  we  first  describe 
a  typical  test  in  some  detail  and  then  describe  results  of  other  one¬ 
dimensional  simulations.  The  reader  will  appreciate  that  completely  cate¬ 
gorizing  the  behavior  and  misbehavior  of  an  approximate  but  complex  com¬ 
putational  model  is  an  extensive  undertaking.  Thus  not  all  possible  tests 
are  shown  or  in  fact  have  yet  been  performed.  We  hope  the  reader  will 
become  aware  of  some  of  the  inherent  levels  of  inaccuracy  in  this  kind 
of  calculation  from  studying  the  benchmark  case  presented  below.  Final 
solution  accuracy  of  better  than  20-30%  seems  unlikely  and  much  larger 
errors  are  possible  when  dealing  with  borderline  phenomena  such  as  critical 
detonation  and  ignition  energies.  Nevertheless,  we  present  these  calcu¬ 
lations  as  the  most  accurate  of  the  type  yet  produced  for  a  specific  physi¬ 
cal  system. 

The  primary  benchmark  test  is  a  simulation  of  a  bursting  diaphragm 
problem  in  Cartesian  geometry.  The  driver  gas  is  and  the  reactive 
mixture  is  a  stoichiometric  l^-air  mixture  at  300°K  and  1  atm  of  pressure. 
The  system  was  10  cm  long  with  the  diaphragm  at  1  cm.  The  analytic 
Rankine-Hugoniot  solution  to  the  problem  indicates  that  if  there  were  no 
viscous  losses  and  no  chemical  reactions,  an  initial  shock  of  Mach  number 
4.3  would  propagate  into  the  detonable  mixture  at  1.7x10 5  cm/ sec.  The 
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temperature  and  pressure  behind  the  shock  front  should  be  1280°K  and  22  atm, 
respectively. 

Figure  4  shows  calculations  of  the  maximum  temperature  in  the  system 
as  a  function  of  time.  The  detailed  calculation  was  carried  out  for  500 
steps  and  the  induction  parameter  model  for  3000  steps.  The  adaptive 
gridding  algorithm  used  110  cells  with  a  minimum  cell  size  of  0.006  cm. 

Also  shown  on  Fig.  4  are  comparisons  of  the  position  of  the  shock  and 
detonation  front  as  a  function  of  time.  Both  of  the  models  show  an  in¬ 
crease  in  velocity  at  about  the  time  that  the  chemical  energy  release 
becomes  significant  and  the  differences  in  these  globally  determined 
integral  quantities  is  slight. 

Figure  5  shows  the  temperature  as  a  function  of  position  at  steps 
300  and  400.  Energy  release  begins  in  the  shocked  material  which  has  been 
heated  the  longest.  Thus  a  deflagration  front  is  created  which  moves 
across  shock -heated  gas.  The  rate  of  propagation  of  the  deflagration  is 
controlled  primarily  by  the  induction  time  of  the  material.  Eventually 
an  equilibrium  is  reached  in  which  the  energy  release  occurs  at  some  short 
distance  behind  the  shock  front.  It  can  be  seen  from  Fig.  5  that  the 
detailed  model  reaches  this  distance  slightly  sooner  than  the  induction 
parameter  model.  By  step  400,  they  have  both  reached  this  equilibrium. 

Several  sources  of  potential  inaccuracy  in  the  global  induction 

parameter  model  have  been  discussed  earlier.  The  largest  of  these  we 

t 

believe  to  be  the  intrinsic  assumption  that  /  dt/x  «  1  implies  the  onset 

o 

energy  release.  This  approximation,  however,  is  an  improvement  to  the 
single  Arrhenius-type  reaction  approximation  because  it  is  constructed  to 
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be  exactly  correct  whenever  the  temperature  is  constant  throughout  the  in¬ 
duction  period*  As  new  results  and  hence  improved  calibrations  become 
available,  improved  analytic  fits  for  generalized  forms  of  the 
model  will  be  made* 

The  calibrated  induction  parameter  model  for  H^-air  was  next  used  to 
investigate  the  minimum  diameter  and  the  associated  minimum  detonation  energy 
for  spherical  geometry.  The  calculations  simulated  a  spherical  bursting 
bubble  problem  with  an  driver  gas  at  1100®K  and  about  700  atm  of  pres¬ 
sure.  Tests  were  performed  in  which  the  volume  containing  the  driver  gas 
was  increased  until  a  minimum  ignition  radius  for  detonation  was  determined. 
This  corresponds  roughly  to  adding  energy  to  the  system  until  a  minimum 
detonation  energy  is  determined.  As  shown  in  Fig.  6,  the  minimum  radius 
for  detonation  is  between  4.9  and  5.0  cm  corresponding  to  energies  between 
1.0  X  106  and  1.1X106  J.  Currently  a  systematic  study  using  this  method 
for  determining  minimum  energies  of  ignition  is  being  carried  out  for 
air  and  CH^-air  mixtures  in  which  the  composition  of  the  driver  gas  and 
the  delivery  of  the  initial  added  energy  are  varied. 

Several  calculations  and  tests  have  been  performed  using  the  stoi¬ 
chiometric  CH^-air  model  described  in  the  previous  section.  Because  the 
induction  times  for  methane  systems  are  always  longer  than  for  systems 
at  equal  pressures,  temperatures  and  stoichiometries,  the  temperature  pla¬ 
teau  between  the  shock  and  burn  front  is  correspondingly  longer.  Results 
of  two-dimensional  simulations  which  use  this  model  are  presented  below. 
Detailed  calibrations  such  as  described  for  the  system  are  currently 
being  set  up  for  methane  and  other  fuels  for  which  the  reaction  kinetics 
are  relatively  well  known. 
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Two-Dimensional  Detonation  Simulations 


Two  multi-dimensional  calculations  of  propagating  detonations  have 
been  performed  using  the  FAST2D  computer  code*13*14  Both  used  the  induc¬ 
tion  parameter  model  calibrated  above  for  stoichiometric  hydrogen-air  and 
methane-air  mixtures.  Realistic  detonations  seem  to  settle  into  dyna¬ 
mically  stable  criss-crossing  cellular  patterns. 24 » 25 »26 *27  The  essential 
instability  of  a  planar  or  near-planar  detonation  to  transverse  pertur¬ 
bations  has  been  previously  demonstrated  by  Flckett  and  Wood1  and  by 
Mader.2  Here  we  are  concerned  with  the  nonlinear  limit  cycle  of  these 
transverse  perturbations  in  the  given  gas  mixtures.3 

The  two-dimensional  calculations  were  both  performed  on  a  100  cell 
by  30  cell  Cartesian  grid.  The  system  and  cell  sizes  for  the  methane-air 
calculation  (Ax  ■  0.15  -  1.5  cm, Ay  *  0.15  cm)  were  chosen  to  be  ten  times 
those  used  for  the  hydrogen-air  calculation.  The  longer  induction  times 
for  methane  lead  to  correspondingly  longer  reaction  zones  and  larger 
stable  detonation  cell  sizes.28  To  permit  spatially  extensive  calculations 
a  region  of  50  fine  zones  centered  on  the  shock  and  driving  deflagration 
are  adaptively  embedded  in  a  grid  of  5C  coarse  cells*  The  total  system 
length  was  82.5  cm.  Solid  wall,  reflecting,  free  slip  boundary  conditions 
were  applied  at  all  four  boundaries:  x*0  and  82.5  cm,  y*0  and  4.5  cm. 

The  result  is  the  two-dimensional  planar  geometry  shown  schematically  in 
Figs.  7  and  8. 

The  system  was  initialized  by  assuming  that  an  oblique  Mach  5  shock 
existed  in  a  small  region  at  the  left  hand  end  of  the  tube.  The  shock 
front  was  oriented  at  60°  from  the  horizontal  x  axis,  as  shown  in  Fig.  8. 


L3 


In  the  shaded  region  from  x»0  to  x*1.5  cm,  the  induction  parameter,  which 
satisfies  the  equation 


vx  3x  y  3y 


1 

t(T,P)  * 


(2) 


was  linearly  ramped  from  1.2  to  0  to  provide  a  smooth  initiation  of  energy 
release.  I(x,y)  was  set  to  zero  throughout  the  remainder  of  the  computa¬ 
tional  mesh.  For  the  hydrogen-air  calculation  a  Mach  4  shock  at  85°  was 
chosen  so  that  we  may  study  the  growth  of  transverse  structures  which 
start  from  the  relatively  small  perturbation  shown  on  the  left  of  Fig.  7A. 
In  both  figures  the  primary  shock  location  is  plotted  at  several  different 
times.  The  number  of  computational  timesteps  is  also  indicated.  The 
narrow  lines  behind  the  primary  shock  at  several  of  the  times  show  the 
instantaneous  location  of  dynamic  pressure  maxima  (secondary  shocks)  in 
the  heated,  shocked  flow. 

In  Fig.  7  the  hydrogen-air  system  displays  growth  of  transverse  per¬ 
turbations  from  a  rather  small  initial  amplitude  which  reach  a  nonlinear 
limit  cycle  by  9.5  usee.  The  temperature  and  pressure  increase  arising 
from  shock  reflection  at  the  upper  and  lower  walls  leads  to  more  rapid 
local  chemical  energy  release.  This  in  turn  enhances  the  reaction  rate 
behind  a  cylindrically  expanding  Mach  stem  which  intersects  the  weakening 
incident  shock.  When  the  Mach  stem  reaches  the  other  wall,  another  reflec¬ 
tion  is  generated  and  the  roles  of  incident  shock  and  Mach  stem  are  inter¬ 
changed  . 

When  the  system  size  is  larger,  as  in  Fig.  7B,  a  more  complex 


pattern  consisting  of  two  detonation  cells  develops.  This  indicates 
that  the  .45  cm  half-width  of  the  upper  calculation  must  be  close  to  the 
stable  free  space  half-cell  size*  Taki  and  Fujiwara3  mention  analogous 
results  for  a  simplified  system. 

Figure  8  shows  a  similar  calculation  for  the  methane-air  mixture 
where  the  initial  shock  is  given  a  larger  angle  from  the  horizontal  to 
speed  attainment  of  a  stable  limit  cycle.  Again  the  system  settles  into 
the  single  mode  fundamental,  suggesting  that  the  4.5  cm  half-width  is 
again  close  to  the  free  space  value.  In  fact,  this  value  agrees  qualita¬ 
tively  with  Urtiew's  estimate.28  The  axial  replication  length  of  this 
pattern  is  about  18  cm,  but  this  clearly  depends  on  the  specified  system 
height,  4.5  cm. 

Taki  and  Fujiwara3  show  that  the  average  propagation  rate  of  the 
detonation,  even  in  its  multidimensional  form,  is  close  to  the  Chapman- 
Jouguet  velocity.  This  is  not  surprising,  since  only  global  properties 
of  the  undisturbed  gas  are  involved.  The  peak  pressure  in  these  two- 
dimensional  calculations  fluctuates  from  a  maximum  at  shock  reflection  to 
a  value  more  than  a  factor  of  two  smaller  just  before  the  next  shock  re¬ 
flection.  This  maximum  occurs  at  the  triple  point  and  decays  continually 
until  the  next  shock  reflection  occurs.  The  characteristic  induction 
distance  at  the  rear  of  the  incident  shock  just  before  reflection  is  much 
larger  than  the  induction  distance  behind  either  of  the  shocks  when  the 
Mach  stem  is  short.  Thus  this  maximum  distance  must  be  intimately  related 
to  the  preferred  structural  size. 
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VI.  Conclusion 


In  this  paper  we  have  described  an  induction  parameter  model  for  use 
in  detonation  and  other  reactive  shock  problems.  The  calibration  of  this 
simplified  model  has  been  discussed  and  results  have  been  presented  for 
one-  and  two-dimensional  problems. 

The  results  described  show  the  importance  of  developing  calibrated 
phenomenological  models  for  use  in  complex  reactive  flow  problems.  Table  I 
clearly  shows  what  may  be  gained  in  computational  time.  However ,  since  the 
programs  required  to  integrate  the  ordinary  differential  equations  descri¬ 
bing  the  chemical  reactions  and  all  of  the  information  about  reaction 
rates  and  intermediate  species  no  longer  must  be  stored ,  a  greatly  increased 
system  size  also  can  be  considered.  This  is  helped  by  the  use  of  the  adap¬ 
tive  gridding  techniques  mentioned  above. 

Thus  we  feel  that  reasonably  reliable  calculations  may  now  be  per¬ 
formed  at  relatively  little  expense  for  one-dimensional  systems  that  may 
be  meters  long.  The  other  gain  for  multi-dimensional  systems  is  the  ability  to 
represent  even  more  complicated  chemistries  inexpensively  by  the  same 
procedure.  Work  is  continuing  on  the  preferred  cell  size  problem  and 
the  problem  of  understanding  cell  regeneration  in  expanding  detonations. 
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Table  I 


Computational  Information  on  Benchmark  Tests 


General  Description 

Detailed  Simulation 

Induction  Simulation 

Number  of  cells 

110 

110 

Distance  simulated 

2  cm 

2  cm 

Minimum  cell  size 

0.006  cm 

0.006  cm 

Maximum  cell  size 

0.1  cm 

0.1  cm 

Computer  Storage 
Requirements 

80  K  bytesa 

b 

Timing  Breakdown 

Number  of  timesteps 

2000 

2000 

Typical  Timestep 

7. 7x10“ 9  sec 

7.7X10'9  sec 

Convective  Transport 

20  sec 

14  sec 

Chemistry 

130  sec 

6  sec 

Total  Time 

150  sec 

20  sec 

a.  This  includes  the  full  software  to  do  the  diffusive  transport  calcu¬ 
lations  which  are  not  included  in  the  work  described  in  this  paper. 

b.  The  induction  parameter  model  was  set  up  as  an  option  in  the  full 
simulation.  Eliminating  the  chemical  integration  routines  would 
greatly  reduce  core  requirements. 
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Table  II 


Induction  Time  Fits 

„  .  a* 

=  xJTi  T*  +  t-t*  *  y  *  ^  p(at®) 

Z  =  y  +  1.25 
H2:02:N2/2:1:4 

■  -  14.64  +  1.30Z2/(6.09+Z2) 

A*  =  4100  -  3850Z2 / ( 7 . 33+Z2) 

T*  *  300  +  1350Z  / (18.75  +  Z  ) 

CH4:02:N2/1:2:8 

i/l  t*  «  -  216  +  0.967y  +  0.104y2  -  0.08y3 

A*  «  2.49X104  -  5.01xl03y  -  4.36xl02y2  +  3.96xl02y3 
T*  *  -  184  +  29 ly  +  23. ly2  -  21. ly3 
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POSITION  (cm) 


Fig.  1  —  Calculations  of  temperature  as  a  function  of  position  in  a 
laser  initiated  detonation.  Three  sizes  of  computational  cells  were 
used  giving  progressively  better  resolution.  The  vertical  line  indicates 
the  correct  position  of  the  bumfront. 
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Fig.  3  —  Calculations  of  temperature  as  a  function  of  time  for  various  initial  temperatures 
at  5  atm  in  a  stoichiometric  CH^air  mixture.  The  arrows  indicate  the  time  of  maximum 
oxygen  atom  concentration. 
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